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On the Numerical Solution of 
Parabolic Partial Differential Equations 1 

Gertrude Blanch 

The numerical results presented here relate to a two-dimensional parabolic partial 
differential equation containing a nonlinear term. Denoting the independent variables by 
t and x, a lattice is introduced, with intervals k and h in the /-and x-directions, respectively. 
Much attention has been devoted recently to the study of the conditions on tne mesh ratio, 
k/h 2 , under which an approximation by a difference equation converges to the solution of the 
differential equation for sufficiently small h. Some known results are summarized in sections 
1 and 2, and three approximation formulas are given, one of order two, and two of order four. 
The feasibility of using approximation formulas of order higher than the differential equation 
is studied in later sections. The primary objective of this paper is to seek the most economical 
mesh ratio for a given approximation formula, that is, of all mesh ratios that will lead to a 
preassigned upper bound of error in the approximation, to choose that mesh ratio that will 
lead to the least amount of work. It is shown in section 3 that the largest admissible mesh 
ratio is not necessarily the most economical one and that a great deal depends on the form 
of the differential system and the boundary conditions. 

In seel ion I -a generalization is given of the method of Hartree and Womersley (1937) 
for improving a solution from two difference approximations. The method is shown to be 
very effective for suitable boundary conditions. Five numerical examples are presented and 
analyzed in section 5. An appendix, with detailed derivations of the formulas used, is given 
for the benefit of those who may want to apply the formulas to specific studies. 



1. Definitions; Basic formulas 

Once the existence of a unique sol tit ion to a dif- 
ferential (^illation has been established, and an 
approximating function has been found that con- 
verges to the solution under suitable conditions, there 
remains the problem of providing an effective numeri- 
cal treatment of the approximation. Our study 
concerns itself with one phase of this problem, for 
the case when the approximating function is ex- 
pressed by a difference equation. We shall further 
limit the discussion to a specific type of differential 
equation, namely, 

with given initial and boundary conditions. Let us 
introduce a lattice covering the region, at intervals 
h in the ^-direction and k in the ^-direction. Let 
\=k/h 2 be the mesh ratio. It is a trivial restriction 
to assume that x a =sh, where s is an integer. 

Notation. For the sake of brevity, we shall write 
u(x,t) = u m ,n, if CM) is a lattice point. 
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where g is any function under consideration. Fur- 
thermore, following the usual convention for even 



1 The preparation of this paper was sponsored (in part) by the Office of Naval 
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central differences, we define, 

w=p / 2r) \ 



m-\-w,n» 



(2) 



The point in the x,£-plane with coordinates x=mh, 
t=nk will be denoted by (m,7i). 

Let u be regular at (m,n). Then for k sufficiently 
small there is a Taylor series in t around (m,n) 



"'TO, w + 1 — ^TO, n I / i 



k"d»u m , n 



^ip! i>t» 



(3) 



If in (3) terms involving p>2 are dropped, and if 
du/dt is replaced by the right-hand side of (1), with 
d 2 u/dx 2 approximated by central differences, we get 
the well-known approximation to u: 

^,B|l = ^ffl,»T^ m ,n \~kj m%n 

= (l — 2\)v m ,n + Hv m -l, n +V m + l,n) + kf m ,n- (4) 

Similarly, if terms through p=2 are retained, we 
obtain, 

2Vn + l = ^TO,tt+ \5 2 + (-^- — — J 5 4 V m!n + k(p m>n 

(X 2 X \ 
-9- — J^ ) (v m -2,n+Vm+2,n) + k(p m ,n, (5) 
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where 



i — l^jm, n \ 



w( 



0- J m, n i ^ J n 



dt 



dx' 



')■ 



Formulas (4) and (5) need modification at the bound- 
ary to satisfy given initial and boundary conditions 
to a required accuracy. This modification can, in 
general, be made. In the present study we shall 
assume that all required derivatives exist and are 
continuous. For parabolic equations this is not a 
serious restriction, for let us assume that we have 
generated values of v for a given t. From (4) and 
(5) it is clear that these serve as boundary values in a 
subdomain for generating the set of values for the 
next t in the lattice. In imposing the continuity 
restrictions, we therefore merely imply that regions 
close to the given boundary, which may have "cor- 
ners" or other discontinuities, will be treated sepa- 
rately. This in fact must usually be done in practice, 
either by choosing a lattice that is much finer than 
that required over the major portion of the region or 
by special approximations that are appropriate for 
the particular problem. Our main concern here is 
with the choice of the mesh ratio, X, for the major 
portion of the domain where the function is presumed 
to be regular. The manner in which an error (or 
variation) in the boundary conditions is propagated 
over the rest of the domain must, of course, be ex- 
amined; but this problem is no different at the 
boundary of the given domain than at any of the 
subdomains (that is, at the successive values of t in 
the lattice). This problem will not be considered 
here. 

Truncation terms. Let <r mtn =u mfn — v m , n - It can 
be verified that corresponding to (4), 

where 

T 2 , m , n =h* [(|-j^) ^gr + ^] + W- (7a) 
Another form for T 2 , m<n is given below: 

T 2 . m .»=^ 2 [|^gF-^^^]+0(P). (7b) 
Corresponding to (5), 

?V?7 Gm,n\~Kl-\ % m,n<, (8) 



where two expressions for jT 4fW(W are given below. 
T 4tm , n =h^M(\) d 



dx* 



-0 2 ]+O(/i 5 ), (9a) 



(9b) 
where 



M ^=J-T2+W 



(10) 



In the above di and 6 2 depend on/ and its derivatives; 
if /=0, then 6 1 = e 2 =0. Further, bf/du and blj>/du 
imply evaluation of these functions at (m,ri) corre- 
sponding to a value of u intermediate between u mtn 
and v m , n . Again, (6), (7), (8), and (9) need modifica- 
tion in the immediate vicinity of the boundary. 

Definition. T Ttm>n will be said to be of order r if 
it contains h r as a factor, but not h r+l . 



2. Stability Considerations; Bounds for the 
Errors in the Approximations 

Definition. The solution v(x,t) of an approxi- 
mating difference equation will be defined as stable 
if it is bounded for all finite t, independently of h. 

Consider the special case when j(x,t,u) = 0, and 
u(x,0) =A{x) is denned and bounded for — oo <Cx<^ °° . 
It will be convenient to refer to the differential eq (1) 
under these special conditions as the basic homo- 
geneous equation. Let this basic homogeneous 
equation be approximated by a difference equation of 
order 2p, so that we can write, 



V m ,n+1 — Z—ib w b 
w=Q 



L m. n — 



2—i®"u$m-\-w,ni 



(ID 



where the coefficients b w and a w are suitable constants. 
It is easy to show that 2& w =l; in general, the 
coefficients a w will be functions of X. Let q be an 
upper bound of |vl(x)|. If it is possible to choose X 
so that all the coefficients a w are nonnegative, we 
shall have, from (11), 



<2JaJfl 



m~\-w, 0, 



<QjZa w =q, 



and by induction on n we can establish that | v Mi n \ < q 
for all m } n; hence v(x,t) is stable. The condition 
that all a w be nonnegative is therefore a sufficient, 
though not a necessary condition 2 for the stability 
of v(x,t) in the case of the basic, homogeneous 

* The above definition of stability and the observation that (1 1) is stable when all 
the coefficients a w are nonnegative were mentioned by F. John in seminar talks 
at the Institute for Numerical Analysis. 
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(Hj nation. We shall refer to the range of X for which 
the basic, homogeneous equation is stable as the 
"admissible" range of A. In what follows X will 
always be chosen to lie within the admissible range; 
this may not guarantee thai the solution v(x,t) will 
be stable for arbitrary boundary conditions and func- 
tions f(x,t,u). However, the choice of X within this 
range usually simplifies the error analysis for any 
set of boundary conditions, even when (1) is not 
homogeneous. 

Turning now to the first approximation formula 
defined by (4) , it is clear that the admissible range of 
X is 0<X<^. The corresponding equations for <j m%n 
are given in (6). Let r be an upper bound of | T 2 , m ,n\ 
and let co be an upper bound of \df/c>u\, this upper 
bound to be independent of h. If o- m , =0, eq (6) 
yields, for a wide range of initial and boundary 
conditions, \a m j\<CkT. Since X is in the admissible 
range, we have further 

Wm,n+l\<Wm,n\(l+fC(X>)\+kT. 

Now (1+ kw)^^'", where t=Nk; hence at time /, 



Wm.Nl^tre* 



(12) 



Boundary conditions may impose modification of 
(12). If the conditions are such that no negative 
power of h is added as a factor to the right-hand side 
of the inequality (12), then it can be shown from the 
form of T 2 , m ,n that an upper bound of r can be found 
that has h 2 as a factor. It follows then from (12) 
that v m , n -^u min as h-^Q. 

Let us now consider the second approxima- 
tion formula, defined by (5). The coefficients of 
v m+Wt n[±w=Q,l } 2] in (5) will all be positive if O^A^f. 
The approximation v(x,t) to the basic homogeneous 
equation will therefore be stable for this range of X, 
and by the same analysis as before, we can show 
that v(x,t) approaches u(x,t) as h approaches zero, for 
a wide range of boundary conditions. 



3. Criteria for the Choice of a Suitable 
Mesh Ratio 

The error a m>n is a function of h,\, and of the bound- 
ary conditions associated with the differential equa- 
tion. Given an upper bound of error that can be 
tolerated in the solution, the problem involves choos- 
ing h and X (the latter within the admissible range) 
so as to meet requirements with the least amount of 
work. We shall assume that for a given approxima- 
tion scheme, the work is proportional to the number 
of lattice points at which u m , n must be evaluated. 
This is not strictly true. For let us define a profile 
as a set of values of v m<n for a fixed n, and all m<s. 
If, for example, successive profiles are generated from 
preceding ones on an IBM machine such as the card 
programmed calculator, then merging operations 
may be required at the end of a profile which may 
consume some time. Thus a grid of 10 points in the 



^-direction and 100 points in the ^-direction may take 
more time to generate than a grid of 20 points in (he 
x-direction and 50 points in the /-direction. Never- 
theless, the assumption that (he work is proportional 
to the number of lattice points is close enough (o 
reality to be useful. Of course, the complexity of the 
programming must be considered; thus an approxi- 
mation formula of order four may take more machine 
cycles (hence more time) than one of order two. 
However, for the same approximation scheme, tin; 
choice of X does not change the amount of work 
radically. Let X=sh be the range of x and ti = Nk 
the range of t. The number of lattice points in the 
region is Ns=(Xt l /hk) = Xt 1 /\h?. As Xt { is fixed, 
the work required for a given approximation scheme 
is therefore inversely proportional to \h. s 

If the exact solution for <r m<n were known, it would 
be theoretically possible to study the magnitude of 
the error for various choices of X and h corresponding 
to a given approximation scheme. The precise 
solution o mM is not easy to find. However, from (12), 
it is clear that an upper bound which can be approxi- 
mated has |r r> m,w| as a factor. We shall, therefore, 
aim to choose X and // in such a manner as to make 
\ r J\, m<n \ small. Moreover, for both approximation 
formulas (4) and (5), the successive terms of h 2 T r m n 
involve h r + p d r + p u/dx r+p and h r + p d r+p f/dx r + p . We shall 
require (hat for all choices of A, (he interval h be 
sufficiently small so as to satisfv 



\8 r+ vu\>c\5 r + p + l u\, />>(), 0<c<l, 



(13) 



almost everywhere. From the known relations con- 
necting derivatives with differences, it is clear that 
(13) implies that successive terms of T r%m # will be 
numerically smaller (ban preceding ones. The 
phrase "almost everywhere" for (he condition (13) 
needs explanation. It may happen that in a region 
where 8 T + P u changes sign, a few entries of 8 r + p u may 
be numerically smaller than corresponding entries 
in the higher differences. Such a case may also arise 
near critical points. In particular the condition (13) 
shall be satisfied by the initially given values and 
-f(x,0,u). In practice one often requires that cbe^or 
i. For one powerful check on the accuracy of com- 
puted values is obtained from the pattern of succes- 
sive differences of the entries. Hence, even if the 
criterion (13) were unnecessary from the viewpoint 
of estimating an upper bound of error in the solution, 
it would still be a desirable condition to impose, in 
order to insure that the computed values difference 
with reasonable ease. We shall further require that 
the term of T Ttnitn involving the lowest power of h 
shall approximate the magnitude of T r<m>n to within a 
factor of two. The fact that a restriction is thereby 
imposed on h must be clearly kept in mind. In 
general a value of h small enough to satisfy (13) will 
not necessarily make \T Ttrriin \ small enough to meet 
requirements for a given upper bound of error in 



3 It is no serious restriction to consider the range t\2& an integral multiple of h. 
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the solution. In the instance when it does not, the 
problem posed is to choose X and h judiciously, so 
as to bring the truncating error within permissible 
bounds with the least amount of work. But if the 
required accuracy is rather low, it may well happen 
that after h is chosen small enough to satisfy (13), 
the error measured by T r>W)W may already be small 
enough to meet all requirements. In that case the 
largest X within the admissible range will, of course, 
lead to the least amount of work. 

With these observations we shall now attempt to 
study the dependance of the truncation term on X 
and on h. Ideally, it would be desirable to classify 
differential equations into several types, according to 
the value of X which is appropriate for equations 
belonging to the type in question. A complete 
classification is difficult to set down, but an attempt 
in this direction is made by considering two types: 

Type I. This type is characterized by the follow- 
ing conditions: 



drdx 2 



dx* U * 



(14) 



When (14) holds, it follows that 



dt 2 ~~ 



d 4 u. 



d 3 u 



c^u 
dx 6 ' 



It can be shown that for systems belonging to type I, 
the terms of T Ttm>n) r<4, do not involve f(x,t,u) or its 
derivatives. The basic homogenous equation belongs 
to this type. 

Type II. This type is characterized by the con- 
dition that successive derivatives of u(x,t) with re- 
spect to t are known to be much smaller than corre- 
sponding derivatives with respect to x (usually of 
twice the order) with which they are associated in 
truncation terms such as (7b) and (9b). It is quite 
easy to find systems belonging to this type; the 
numerical illustration given in section 5 belongs to 
type II. 

Consider differential systems of type I, and let us 
start with the case when formula (4) is used. The 
corresponding truncation term, T-, my n, is given by 
(7a) . It has been observed by Milne, as well as by 
Thomas and others, that if X is chosen as -g-, then the 
leading term of T 2>m , n vanishes, so that the order of 
T 2 , m ,n becomes four. With \=\ and h small enough 
to satisfy the conditions for an upper bound of 
I T 2 , m ,n\ — enough to meet the requirements for a given 
upper bound of error in the solution — a certain 
amount of work will be done, which we shall measure 
in units of Z=l/\h 3 , as it has been shown that for 
a rectangular lattice the work is proportional to Z, 
approximately. If another X and h were chosen, 
then we shall ask whether, for the same amount of 
work 7 and using the same formula (4), greater ac- 
curacy can be obtained in the solution, assuming 
that accuracy to be measured by the magnitude of 



the truncation term. In other words, we shall com- 
pare various choices of X and h, for a constant Z. It 
can be readily shown that for equations belonging to 
type I, no other choice of h and X, keeping Z constant, 
will be as good. In this sense it is correct to state 
that X=-| is the best choice for formula (4). How- 
ever, one important point has been overlooked. We 
have seen that h is not completely free, for h must be 
small enough to meet the minimum conditions im- 
posed by (13). It may happen that after h has been 
taken small enough to insure that (13) is satisfied, 
the conditions for an upper bound of error in the 
solution can be satisfied with some range of X>^. 
In that case, we would certainly get a more accurate 
solution by choosing X=-^, but that would be more 
accuracy that required, and work could be saved by 
choosing a larger X. We conclude that even for 
equations of type I, the choice of X=^ will be best 
only if a relatively high accuracy is required in the 
solution. Moreover, in practice the differences of 
u(x,t) can usually be judged only from the initially 
given profile, at the time when X and h are chosen, 
so that a safe interval h (rather than the maximum 
permissible for the given profile) is often chosen, and 
it may happen that the value of h considered neces- 
sary may, as stated, bring the error within the 
tolerance limits for all choices of X. 

Still considering systems that belong to type I, let 
us now examine approximations of order four. The 
expression for the truncation term T 4 , mf/i is given in 
(9a) and the term of order four is 



h*At(x,t,\)-- 



<- 



A . 1 \OU mn i±i\/f/\\ O l^m, n 



It can be readily verified that M(X) is positive for all 
choices of X, hence the leading term of T 4>mi7l cannot 
be eliminated completely, as in the case of the 
simpler approximation of order two. However, we 
may seek that value of X which will make the unit of 
work, Z=l/\h 3 , a minimum, subject to a given 
permissible upper bound of \T it7rltn \. By hypothesis, 
h is small enough so that h*A±(x,t,\) approximates 
the magnitude of | T itm>n \ ; hence if C is the permissible 
tolerance of |T 4>OT , n |, we wish to satisfy 



A 4 |^(s,*,X)|=A 4 Af(X) 



cfu 



dz ( 



<C } 



(15) 



Since \d 6 u/dx 6 \ is independent of the choice of 
parameters h and X, the inequality expressed in (15) 
can be satisfied by taking 



A 4 M(X)=<7 X ; Ci<C/ 



dx ( 



(16) 



Thus we seek to determine h and X, which minimize Z, 
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subject to the condition 

h 4 M(\)=Ci. 



(17 



From (17), h=C 1 /[M(\)]K Substituting this value 
of A into Z and differentiating the resulting expression 
of Z with respect to X, we obtain the following condi- 
tion for a minimum of Z: 



or 



F=3XM'(X)-4M(A)=0, 
60X 2 +15X— 8=0. 



(18) 



The positive root of (18) is \ =0.26095 . . ., and it 
can be verified that F (\) is positive, so that Z is 
indeed a minimum for this value of X. Since X lies 
within the admissible range, it can be used for the 
mesh ratio, in conjunction with a suitable value of h 
which satisfies (18). In practice, X=0.25 will 
normally be used, because an irrational value of X is 
inconvenient. 

It will be instructive to examine the following 
schedule, which gives M(X), h, and the work unit Z 
for a constant value of A 4 M(X), and various values of X 
within the admissible 1 range, in units of the corre- 
sponding quantities when \=i. 



X 


*«-(?-n+55) 


h 


Z=l/(X/> 3 ) 


% 


0. 001851 


1. 565/ii 


0. 782.-1 


}i 


. 0011 11 


1. 778/m 


. 444-! 


Y. 


. 000694 


2. 000/m 


. 25-i 


Mo 


. 001111 


1. 778/n 


. 297:, 


% 


. 001851 


1. 568/i.i 


. 391:, 


% 


. 003299 


1. 354/ii 


. 536;, 


V?. 


. 011111 


h 


Z\ 


% 


. 029630 


0. 783/ii 


1. 565*1 



Compared with a unit of work Z x for the ease X=|, 
only \Z X is required when X=|. The largest admis- 
sible value of X is the poorest of all, from the view- 
point of the amount of work required. Again, a 
word of caution is required. For the same magnitude 
of the error in the leading term of T± tm%T and different 
choices of X, the above schedule shows that with 
X=J, h can be chosen twice as large as that required 
when X=§. If A is made twice as large, the eighth 
difference in the /-direction is multiplied by about 2 8 , 
and the differences must be reexamined to see whether 
this larger value of the eighth difference still satisfies 
the fundamental conditions imposed by (13) namely 
that successive differences in the ^-direction beyond 
the fourth be numerically smaller than preceding 
ones. Moreover, it has been pointed out before 
that the maximum h corresponding to which (13) 
is satisfied may already be such that T^.m.n is within 
the required tolerance for all admissible values of X. 
In that case it is of course best to choose X=f or one 
close to it. 



The following question arises: since the simpler 
approximation (4) has a truncating error of order 
four when X=4, and T Ajn<n is also of order four, is 
there any gain in using the approximation formula of 
higher order? The answer to (his question is com- 
plicated by the fad that the time required to gener- 
ate a profile corresponding to the higher approxima- 
tion may be considerably longer than that for the 
simpler approximation. Much will depend not only 
on the computing instrument which will be used, but 
also on the complexity of the boundary conditions. 
If the IBM card programmed calculator is to be 
used, and the boundary conditions are not too com- 
plex, the simpler approximation (4) can perhaps be 
generated in only two-thirds of the time per profile, 
compared with the more elaborate fourth order ap- 
proximation given in (5) . There are, however, com- 
pensating factors which make the higher approxima- 
tion worth considering. Let us assume we are dealing 
with a case where the "best" values of X are used in 
the approximation of order two and the one of order 
four. For the same h used in both cases, the number 
of lattice points is inversely proport tonal to X. Hence, 
we shall use only two-thirds the number of lattice 
points when the higher approximation is used. This 
would about compensate for the longer time it may 
take to generate each profile. There is, however, a 
gain in accuracy when the higher approximation is 
used. For the coefficient of h^dH/dx®) in (9a), cor- 
responding to the simpler approximation, is 1/540 
when X=^. On the other hand, the coefficient of the 
corresponding term in the higher approximation with 
X=-^is only 0.000694; or less than three-eighths that 
of the simpler approximation. Furthermore, in 
cases where the maximum permissible h is such that 
Ti, m , n is already within the required tolerance limit 
for all values of X, we may be able to take X=f , or 
close to it, when the higher approximation is used. 
Such a choice of X woidd cut the number of lattice 
points to \ that required for the simpler approxima- 
tion, with X=-^, and that might more than compen- 
sate for the greater difficulty in generating v(x,t) by 
the higher approximation. 

Let us now consider differential systems that be- 
long to type II. Since, by hypothesis, the deriva- 
tives in the ^-direction are negligible compared with 
those in the x-direction, (7b) shows that the leading 
term in T 2 , m , n is practically independent of X, and 
this term cannot be eliminated by any choice of X. 
Hence, it is reasonable to take the maximum admis- 
sible X, or one close to it. Similarly, the leading 
term of T 4tMin is complex in structure, and there is 
no optimum X that stands out as suitable for all 
functions falling under this type. For a given prob- 
lem, it may be possible to compute estimates of the 
various terms that contribute to the truncation error. 
Or if the problem involves a family of parameters, 
in the boundary condition, then available results for 
some members of the family may lead to an optimum 
choice of X and h for the remaining ones. 
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For equations falling under type II the higher 
approximation formula usually has distinct advan- 
tages over the simpler one of order two. For in this 
case the truncation term is necessarily of a lower 
order of magnitude in the higher approximation 
formula than in the simpler one; moreover, a larger 
mesh ratio, namely, §, can be taken. Whenever the 
boundary conditions are such that the coding prob- 
lem is manageable, the higher approximation formula 
is to be recommended. 

It might be well to remark that in some cases the 
solution may be an oscillating function of x but not 
of t. Such functions may fall under type II, and a 
Fourier approximation may actually be better than 
a finite difference approximation. However, there 
is no reason to expect that all, or even a good portion 
of functions belonging to type II, can be most simply 
treated by Fourier approximations. A finite differ- 
ence approximation such as (4) or (5) is often pre- 
ferred, because of its simplicity, to other types of 
approximations (by Fourier series or perhaps "im- 
plicit" difference approximations). Our concern 
here, as stated before, is to study the dependence of 
the solution on the choice of A, after either (4) or (5) 
has been selected as the approximation formula. 

4. Method of Improving the Solution from 
Two Difference Approximations 

The method to be explained below has been 
presented in [1] 4 by Hartree and Womersley, who 
ascribe the idea to L. F. Richardson [5]. In [1] the 
method is applied to a somewhat different type of 
approximation — a mixed difference — differential 
scheme, suitable for computation by a differential 
analyzer. The results will here be extended to 
difference approximations of any other. 

Let us suppose that values of v m>n have been gen- 
erated by an approximation formula, corresponding 
to a true solution u mtn , and let us assume that it is 
possible to write 



0" m, n wm.n v m.n 



--h r C r (x, t, \)+h r ^ l C ri . 1 (x ) t, X)- 



(19) 



where the functions Cj(x,t,\) are independent of h. 
It can be shown that corresponding to a wide range 
of boundary conditions the expression for <j m>n does 
assume the form (19) for both approximation 
formulas (4) and (5). 

Consider now the case when v m , n has been gen- 
erated by the use of an interval h x in the ^-direction 
and a mesh-ratio X. Let these values of v mtn be 
designated by fl m>n (Ai). Now let another computa- 
tion be made, based on the same mesh ratio X, but 
at an interval h 2 in x, where h 2 = ph 1 , 0<^p<^l. These 
values will be designated 5 by v m>n (h 2 ) . It is presumed 
that the same approximation formula will be used in 



4 Figures in brackets indicate the literature references at the end of this paper. 

5 In the subsequent discussions we shall write v(h,2), or v(hi,h2) to indicate 
m.nihz) or Vm, n (hi,h2), when no ambiguity is likely to arise. 



both cases. Clearly, 

v m ,o(hi) =v mt0 (h 2 ) =u mt0 ; v Stn (h 1 )=v Sin (h 2 )=u Smn . 
By hypothesis, we have from (19) 



z—v m , n (h 1 )=h r 1 C r +h r 1 +1 C r - 



\m (^r+2~r • • • 



(20) 



Similarly, remembering that h 2 =ph l and that Cj is 
independent of h, 

n m ,n — v m ,n.(h 2 ) = P r h r 1 C r +p r+1 h r 1 +1 C r+l 

+ P r +% +2 C r+2 + ... . (21) 

Multiply (20) by p r and subtract from (21). This 
gives, after transposing some terms and dividing by 
(1-P r ), 



,n=V m>n (h 1 ,h 2 )- P \ {l r P) A r + 1 C r r+l 
1 — O 



where 



■p ri j~h r+2 Cr + 2+. 



v m ,„{h 1 ,h 2 )= 



Vm,n(h 2 )—p T V m , n (h,) 



l-p r 



(22) 
(23) 



We shall refer to the process defined in (23) as the 
"p T " correction. In [1] it is recommended that p be 
taken as ^. This is a desirable choice for numerical 
work, since every other value in the ^-direction at the 
smaller interval is available at the larger interval. 
Similarly, every fourth value in the ^-direction will 
be available at both intervals. There are many ways 
of applying corrections of this type. One way, for 
example, is to generate four values in the ^-direction 
at the finer interval, then generate the corresponding 
values at the larger interval; apply the p r correction 
to the last profile, and use these corrected values as 
initial data for generating the next profile. Such a 
use of the correction scheme would, of course, require 
interpolating for values of v(hi,h 2 ) for every other 
value of x, since we can correct only for those points 
that are available at both intervals. The coding of 
this method would be complicated. The simplest 
way of using the correction scheme is to actually 
make two completely separate computations for all 
required values of n, and then to apply the correc- 
tion process only to those functional values that are 
required. Often the last profile generated is of most 
interest, and in that case it may be enough to correct 
the values on the last computed profile only. 

The process furnishes a powerful check on the con- 
vergence rate of the approximations, and since the 
work done at the larger interval is one-eighth that 
done at the smaller interval if p= J, the added labor 
is not too costly, when the two computations are 
carried separately. Moreover, the coding is the same 
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for both approximation schemes. For the special 
case when r=2, (22) and (23) reduce to 



^m, n — Vm, n\h\Jl 2 ) 



1 + P " 



•p 2 /i 4 C 4 + 



(24) 



where 



/t t n V m , n (h 2 )—p V m , n (hi) , v 

1 — P 



The "p T " correction has been applied in the numerical 
examples given in section 5 with highly satisfactory 
results. In [1] Hartree and Womerslev give sufficient 
conditions on the nature of the boundary for the 
method to be valid; the method can probably be 
used over a wider class of functions than those 
specified in [1]. 

If hi is sufficiently small, v(h u h 2 ) will always be an 
improvement over u(h 2 ), since the truncation term 
is of lower order of magnitude in v(hi,h 2 ). The 
question arises: how small must h be to insure that 
v(h u h 2 ) shall be an improvement over u(h 2 )l Let 
the truncation term, T r , m ,,» now be identified by 
T T , m> n(h)j when associated with an interval h in the 
or directions, and let 



a(h h h 2 ) = u—v(h u h 2 )) <r(h) 



-Mh)) 



T r . m , n(h u h 2 )= Tt > m ' n{fl2 \ pr ? r > m > n{hl) - (26) 
1 — p 

It is more difficult to get a good estimate of a than 
of T Ttm>n . Just as in section 3, we shall, therefore, 
inquire under what conditions the inequality 

\Tr tm ,nQlM\<\T Ttmt nQl2)\ ' (27) 

will be satisfied. A numerically smaller truncation 
term will usually be associated with a smaller total 
error. For the approximation formulas considered 
here the truncation term is of the form 



T T , m , n (h 2 )=JZp v h v g P {x,t,u). 

Moreover, the condition following eq (13) guarantees 
that 

T r , m ,n(h>2)^P r h r g r (;x,t,u). (28) 

But no similar statement was made about T TtTntn (hi). 
Now from 

-p/(l-p p W + »g r+P (x,t,u) 

^r,Tn,n\hi,fl 2 )= yz zz — 

(1 — P ) 



If we can satisfy 

|A<7r+ P (^M)|<|p C ^r+p-l(«,^)|, 



(29) 



where 



P c =l~P r , 



(30) 



then from (27) and (30) it can be shown by summing 
the absolute values of the terms of r I\, m ,n{h\h 2 ) that 
(27) is satisfied. When /> = ■§■ and r=2, the condition 
(29) implies that p c =\ . If p= J andr=4, then p c =|f, 
or what is equivalent, successive terms of T r<mtn at 
the larger interval are required to be just smaller 
numerically than preceding terms. Under such con- 
ditions v(hi ho) will always be an improvement over 
u(h 2 ). 

In the foregoing, it has been assumed thai r mn can 
be computed exactly by the prescribed formula, and 
that all initial and boundary conditions are exact. 
This can seldom be realized in practice, and there 
will be rounding errors committed at every step of 
the computation, due to carrying a fixed number of 
decimals or sigmficanl figures in the computations. 
The cumulative effect of such errors can perhaps be 
studied statistically, or upper hounds for errors of 
this type [2] can be found. In 1 he numerical examples 
given in section 5, the cumulative round-off error 
was very small, after 100 steps in /. 



5. Numerical Examples 

The problem selected for analysis was the following 
one: 



At 6 = 0, u=(10x 2 +x 3 +0M) H =A(.r); 
at x=l, u=(UM+f) H =B(f). 



(31) 



[SL- 0; fMm 



-(3x+9.5) (10.r+1.5r0 2 



u° 



Tins particular form was chosen because it repre- 
sents a case in which derivatives with respect to / 
are much smaller numerically than corresponding 
derivatives with respect to x, with which they are 
associated in truncation terms such as those of (7b). 
The differential system belongs to type II of section 
3. It is known that the exact solution to the problem 
is 

The choice of a nonlinear form of the differential 
equation was deliberate. It was made in order to 
study the cumulative error in cases where a con- 
siderable number of operations, that involve the 
approximate values of u, have to be performed at 
each step. 

The following five sets of solutions were generated : 
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Example 1. 
Formula: 



Vm, w-f 1 — ^m, n\ ™ ® m,n\ "J m, n-> ^m, o ^\.\X) , 

v s>n =B(t)', sh=l. 



>x_] x= 



^0,71 



Parameters: h=0.05=h 2 ; £=0.00125; X=§. 

Range of t: 0<t<0.125. 

Number of lattice points: 20X100 = 2,000. 

Initial values and the computed values of v ntn for the last profile are given in table 1. 

Table 1. Solution of the difference equation 



V m ,n+l = V m ,n-\-^ 2 V m ,n + kf(x,u)', V- Un = V Un - 



JL 



Att = 0, v{x,0) = (10x 2 +x 3 +0.64)^; Atx=\, v(\,t) = (11.64 + 0* 

., N -(&C + 9.5) , (10x + 1.5^)2 . _ ± , A nR . nA1QK 

f(x,v)= — H -^ — ; \ = | ? ^2 = 0.05; ft = . 00125 



u = ^l0x 2 +x*+t + 0M 









At £ = 0.125 




X 


v(x,0) 






« 












v(x,t) 


u(x,t) 


u (x, t)—v{x,t) 


0.00 


0. 8000000 


0. 8720400 


0. 8746427 


+ . 0026027 


.05 


. 8155520 


. 8864322 


. 8888897 


+. 0024575 


. 10 


. 8608136 


. 9285124 


. 9305912 


+ . 0020788 


. 15 


. 9318664 


. 9950813 


. 9966819 


+ . 0016006 


.20 


1. 0237187 


1. 0819105 


1. 0830512 


+• 0011407 


.25 


1. 1316470 


1. 1848335 


1. 1855905 


+ . 0007570 


.30 


1. 2517987 


1. 3003062 


1. 3007690 


+ . 0004628 


.35 


1. 3812585 


1. 4255420 


1. 4257892 


+ . 0002472 


.40 


1. 5178933 


1. 5584295 


1. 5585249 


+ . 0000954 


.45 


1. 6601581 


1. 6973964 


1. 6973876 


-. 0000088 


.50 


1. 8069311 


1. 8412718 


1. 8411952 


-. 0000766 


.55 


1. 9573898 


1. 9891819 


1. 9890638 


-. 0001181 


.60 


2. 1109240 


2. 1404665 


2. 1403270 


-. 0001395 


.65 


2. 2670741 


2. 2946235 


2. 2944770 


-. 0001465 


.70 


2. 4254896 


2. 4512642 


2. 4511221 


-. 0001421 


.75 


2. 5858993 


2. 6100866 


2. 6099568 


-. 0001298 


.80 


2. 7480902 


2. 7708507 


2. 7707399 


-. 0001108 


.85 


2. 9118937 


2. 9333664 


2. 9332788 


-. 0000876 


.90 


3. 0771740 


3. 0974789 


3. 0974182 


-. 0000607 


.95 


3. 2438210 


3. 2630627 


3. 2630315 


-. 0000312 


1.00 


3. 4117444 


3. 4300145 


3. 4300145 






Example 2. 

Formula: The same as in example 1. 

Parameters: A=0.1=/ii; k=0.005; 

Range of t: The same as in example 1. 

Number of lattice points: 10X25—250. 

Initial values and the computed values of v m , n for the last profile are given in table 2. 



X=*. 
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Table 2. Solution of the difference equation described in table 1 and values of v (hi,h 2 ) 
Parameters: fo=0.1; fc = 0.005, X=h 

v(xSJhJi2)^v(h h h 2 ) = [v(h 2 )-p>-v(h l )]/(l~ P *); P = l 
v(h 2 ) given in table 1. 





v(x,0) 


At £ = 0.125 


X 
















v(x,t) 


u(x,t) 


u(x,t) —v(x,t) 


v(hi,v 2 ) 


u{x,t) — v(h\,h 2 ) 


0.0 
. 1 
.2 
.3 
.4 


0. 8000000 
. 8608136 

1. 0237187 
1. 2517987 
1. 5178933 


0. 8638878 
. 9221349 

1. 0784844 
1. 2989434 
1. 5581507 


0. 8746427 
. 9305912 

1. 0830512 
1. 3007690 
1. 5585249 


+ . 0107549 
+ . 0084563 
+ . 0045668 
+ . 0018256 
+ . 0003742 


0. 8747574 
. 9306382 

1. 0830525 
1. 3007605 
1. 5585224 


-. 0001147 

-. 0000 170 

. 0000013 

+ . 0000085 

+ . 0000025 


.5 

. 6 
. 7 

. 8 
.9 


1. 8069311 

2. 1109240 
2. 4254896 

2. 7480902 

3. 0771740 


1. 8415191 

2. 1408939 
2. 4517066 

2. 7711921 

3. 0976676 


1. 8411952 

2. 1403270 
2. 4511221 

2. 7707399 

3. 0974182 


-. 0003239 

-. 0005669 

. 0005845 

. 0004522 

. 0002494 


1. 8411893 

2. 1403240 
2. 4511107 

2. 7701369 

3. 0974160 


+ . 0000059 
+ . 0000030 
+ . OOOOO. - ) 1 
+ . 0000030 
+ . 0000022 


1.0 


3.4117444 


3. 4300145 


3. 4300145 











Example 3. 

Formula: The same as in example 1. 

Parameters: A=l/14; £=1/1176; \=£. 

Range of t: The same as in example 1. 

Number of lattice points: 14X147=2058. 

Initial values and the computed values of v mtU for the last profile are given in table 3. 

The formula for examples 1 to 3 comes from (4). 

Table 3. Solution of the difference equation described in table 1 
Parameters: fr=l/14, A— 1/1176; A=l/6, x=mh. 



m 


V&O) 


At * = 0.125 


v(x,t) 


u(x,t) 


u(x,t)—v(x,t) 



1 
2 
3 
4 

5 
6 

7 
8 
9 

10 
11 
12 
13 
14 


0. 8000000 
. 8314955 
. 9203245 

1. 0531018 
1.2164087 

1. 4003800 
1. 5985781 

1. 8069311 

2. 0228433 
2. 2446210 

2,4711277 
2. 7015788 

2. 9354176 

3. 1722397 
3.4117444 


0. 8692554 
. 8987549 
. 9824630 

1. 1087660 
1. 2656385 

1. 4438567 
1. 6371338 

1. 8413398 

2. 0537601 
2. 2725815 

2, 4965649 
2, 7248430 

2, 9567941 

3. 1919624 
3. 4300145 


0. 8746428 
. 9035402 

. 9858991 

1. 1108660 
1. 2667479 

1. 4443213 
1. 6372086 

1. 8411953 

2. 0535080 
2. 2722948 

2, 4962917 
2, 7246151 

2, 9566326 

3. 1918811 
3. 4300145 


+ .0053874 
+ . 0047853 
+ . 0034361 

+ . 0021000 
+ . 0011094 

+ . 0004646 
+. 0000748 
. 0001445 
-. 0002521 
-. 0002867 

-. 0002732 

-. 0002279 

. 0001615 

-. 0000813 
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Example 4. 
Formula: 



Vm, »+l — Vm, n I A (p y 2" " ) ® m, n\^Jm, n 



'■Vr, 



ran _ m 3 h s 

Ldx] x=0 - Vmn V Q , n 



m- 



(32) 
1,2. 



It will be convenient to rewrite (32) in the following 
form: 



^flt,«4-l = s i &qV m + q,n\ iCJm 
q=-2 



(33) 



where 



T2 A, 



a_i=ai=-|X; a = 1 — 2.5 X. 



Parameters: A-0.0625; & = 0.0014; X=0.3584. 

Range of,*: 0<£<0.1288. 

Number of lattice points: 16X92 = 1472. 

Formula (32) is a modification of (5). It is ob- 
tained by dropping the terms of (3) involving p>2. 
However, the second derivative with respect to x is 



approximated by differences, including the fourth 
order. The resulting formula is almost as accurate 
as (5) for this problem, since derivatives with respect 
to t are very small numerically. In view of the fact 
that (5) was modified, it is necessary to examine (33) 
to determine the admissible range of X. It is clear 
from (33) that all the coefficients a Q cannot be made 
positive by any choice of X; hence the stability cri- 
teria given earlier do not apply. However, it has 
been show T n in [2] that if there exists a positive num- 
ber M, independent of x and t, such that the coeffi- 
cients a q of (33) satisfy 



i S a« exp (iqy) 

lp— 2 



<exp(-iWy), for|y|<7r, 



(34) 



then the basic homogeneous differential equation, 6 is 
stable. From that result we may then deduce the 
stability of (32) or (33). 



stability arc given in 
Table 4. Solution of the difference equation 

--V m ,n+( A5 2 — j2 ) V m ,n-\-kf mtn 

ra 3 /i 3 |~ d fl m z h 3 

v - m ' n ~ Vm ' n ^^Ldx] x=Q - Vm ' n ~1^7 

& A V 8 -i tn =8 A V s -2,n; Sh=l 



6 In [2] the theorerrTapplies to a more general case. Similar results relating to 
" [4]. 



V m ,n+1~ 



V s ,n=Us 



f(x,v) defined in table 1. 

Parameters: h 2 = 0.0625; /c = 0.0014 ;f\ = 0.3584. 



X 


v(x,0) 


At £ = 0.1288 


v(xj) 


u(x.,t) 


u(x,t) —v(x,t) 


0. 0000 
.0625 
. 1250 
. 1875 
. 2500 

.3125 
. 3750 
.4375 
. 5000 
. 5625 

. 6250 
. 6875 
. 7500 
. 8125 
.8750 

. 9375 

1. 0000 


0. 8000000 
. 8242006 
. 8934221 
. 9990767 

1. 1316470 

1. 2833862 
1. 4487872 
1. 6241314 
1. 8069311 

1. 9955052 

2. 1886961 
2. 3856893 
2. 5858993 
2. 7888957 

2. 9943567 

3. 2020364 
3. 4117444 


0. 8767367 
. 8988926 
. 9627974 

1. 0615954 
1. 1872134 

1. 3326403 
1. 4925905 
1. 6633207 

1. 8422337 

2. 0275254 

2. 2179277 
2. 4125343 
2. 6106854 

2. 8118924 

3. 0157864 

3. 2220847 
3. 4305684 


0. 8768124 
. 8989475 
. 9628100 

1. 0615804 
1. 1871920 

1. 3326211 
1. 4925762 
1. 6633105 

1. 8422269 

2. 0275206 

2. 2179247 
2. 4125324 
2. 6106848 

2. 8118925 

3. 0157871 

3. 2220858 
3. 4305684 


+ . 0000757 
+ . 0000545 
+ . 0000125 
-. 0000150 
-. 0000214 

-. 0000192 
-. 0000143 
-. 0000102 
-. 0000068 
-. 0000048 

— . 0000030 
-. 0000019 
-. 0000006 
-f. 0000001 
+ . 0000007 

+. 0000011 
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Example 5. 



Formula: The same as in example 4. 

Parameters: ^=0.125; fr=0.0056; X=0.3584. 

Range of t: The same as in example 4. Results are given in table 5. 

Table 5. Solution of the difference equation described in table ' t 

Parameters: &i = 0.125; & = 0.0056; /i = 0.3584 

v(x,t,h h h 2 ) = v(h h h 2 ) = [v{h 2 )--f>*(v(hi))]/[l-p*]; P =\ 

v(h 2 ) given in table 4. 



X 


v(x,0) 


At t = 0.1288 


v(x,t) 


u(x,t) 


u(x,t)—v(x,t) 


v(hi,h 2 ) 


u(x,t) —v(h\,h 2 ) 




0. 125 
. 250 
.375 
. 500 

. 625 
. 750 
. 875 

1. OOO 


0. 8000000 

. 8934221 

1. 1316470 
1. 4487872 

1. 8069311 

2. 1886961 
2. 5858993 

2. 9943567 

3. 4117444 


0. 8753628 
. 9623470 

1. 1872934 
1. 4926595 

1. 8422296 

2. 2178929 

2. 6106245 

3. 0157128 
3. 4305684 


0. 8768124 
. 9628099 

1. 1871920 
1. 4925762 

1. 8422269 

2. 2179247 

2. 6106847 

3. 0157869 
3. 4305684 


+. 0014496 

+ . 0004629 

. 0001014 

. 0000833 

-. 0000027 

-f. 0000318 

+ . 0000602 

+. 0000741 




0. 8768282 
. 9628274 

1. 1872081 
1. 4925859 

1. 8422340 

2. 2179300 

2. 6106895 

3. 0157913 


-. 0000158 
. 0000175 
. 0000161 
. 0000097 

-. 0000071 

. 0000053 
. 0000048 
-. 0000044 




In all the above examples, v n 



and u n 



were generated. Fourth differences were generated 
in all the examples, even where the computing for- 
mula did not call for them. These fourth differences 
were jotted down by the operator of the card pro- 
grammed IBM calculator, and any lack of continuity 
in the differences was a warning that the machine 
was not functioning properly. The operations were 
carried out with a board wired to perform eight- 
place multiplication. The calculations were carried 
to the fullest possible accuracy, that is, to seven 
decimals in v mn and to eight decimals in some of the 



subsidiary computations for (i 



m. n-f 1 " 



-Vm.n)' This 



accuracy is far in excess of the truncation error, for 
small values of x. 

A u p 2 " correction was applied to the values of the 
last available profile, as explained in section 4, based 
on the entries in examples 1 and 2. The results are 
given in example 2. Similarly, a p 4 correction was 
applied to the last profile in examples 4 and 5; the 
results are given in example 5. 

5.1. Observations 

(a) The p r Correction Process 

It is to be noted that in spite of the fact that in the 
last profile v 0jTl differs from the true values by 0.0026 
in example 1 and by more than 0.01 in example 2, 
the values of v(h u h 2 ) resulting from the p 2 correction 
are correct to within 0.00011 or better. It can be 



verified from the tabulated entries that 
[u-v{h l ))l[n-i>%)]^K\IK 



It follows that to gain an accuracy comparable to 
that of v(hi,h2) without the a p 2 " correction, an in- 
terval of A = 0.01 would be needed, and hence the 
amount of work would be 125 times that used in 
example 1 . 

The improvement due to the "p 4 " correction in 
examples 4 and 5 is not quite so striking. However, 
even here there is considerable improvement for 
small values of x, and it must be remembered that 
the formula used does include an h 2 term in the trun- 
cation error, which, is not eliminated by the p 4 cor- 
rection, although, the error from this term is somewhat 
lessened. However, for x larger than J, v(h 2 ) in 
example 4 is closer to the true value than v(hi.h 2 ). 
An explanation for this may come from the following 
considerations: Assume 



«-»W=Sp p+2 1 ,+2 ^+2. (35) 

p = 

then it follows that 

_P&=£l h *o t -.... (36) 
1 — p 
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It is to be noticed that the terms involving h 2 and A 3 
are numerically smaller in (36) than in (35). How- 
ever, the term involving h 5 is somewhat larger, if 
P=i, and all subsequent terms will be larger. In 
fact, it can be readily verified that if a "p r " correc- 
tion is applied, then compared with a term involving 
p p+r h p+r C p+r in (35), there is the corresponding term 
[- P r (l- P p )/(l-p r )h p+r C p+r in [u-v(h u h 2 )], which is, 
of course, numerically larger; but usually the leading 
term after the elimination is smaller than the term 
eliminated, namely, h r p T C r . The "p 4 " approximation, 
as applied in examples 4 and 5, is rather unusual in 
that a very small term of order two is left. If all the 
terms of (35) are small numerically, it may happen 
that the combination of the leading terms m\u— 
^(hi,h 2 )\ may be somewhat larger than in \u—v(h 2 )\. 
In such cases, however, the difference between v(h h h 2 ) 
and v(h 2 ) will itself be small; hence, although it may 
not be known which is the better answer, it is to be 
expected that the order of magnitude of the error in 
v(hi,h 2 ) is not greater than \v(h 2 ) — v(h h h 2 )\. 

(b) The Rounding Error 

The last place of u(x,t) is not guaranteed, hence we 
can judge the rounding error only to the extent that 
the sixth decimal place is affected. It is to be noticed 
that in example 4, which seems to be the most ac- 
curate, the difference u m , n — v m>n is systematic as far 
as sign is concerned. This is evidence of the fact 
that no large rounding-off error accumulated, after 
100 steps in t, even though a nonlinear differential 
equation was used, and a considerable number of 
arithmetic operations were performed at each step. 



(c) The Effect of Varying X 

Let us compare the error pattern in example 1, 
where X=|, with that of example 3, where \=f. 
In spite of the fact that somewhat more work was 
performed in generating values in example 3, the 
results are not as good — although the error in both 
examples is in the same decimal place; but the error 
in example 3 is about twice as large. This, in fact, is 
precisely what was to be expected. For as observed 
in section 3, in cases coming under type II, where 
derivatives with respect to t are relatively small, the 
error is not appreciably affected by X; hence, it is 
almost proportional to h 2 . The ratio of the two values 
of h 2 is (20/14) 2 ^2.04. The numerical illustration 
verifies the observation that X=| is not necessarily 
the best mesh ratio for the formula given in (4). 
The boundary conditions and the form of f(z,t,u) 
determine the type that the differential system be- 
longs to, and it is only after the problem is studied 
from this viewpoint that a suitable choice of X can 
be made. 

In connection with example 3 corresponding to 
X=^, it was not desirable to generate values at double 
the interval in the ^-direction, for purposes of apply- 
ing a " p 2 " correction. For if the formula given in 
(4) were used to compute 0_i,i, we would get a nega- 
tive value of (y Q) i—v 0t0 ). Even if the true answer 
were not known, a knowledge of the more accurate 
values of v m , n from the smaller value of h would warn 
us that a negative value of (v 0d — v Qt0 ) is incorrect. 
It will be instructive to write down the differences 
of v m ,o (the known initial profile) at double the 
interval used in example 3, and to use a value y_ 1(0 
computed from formula (4). The values of v mi0 and 
some of the differences are given below: 







h 


=4; t=0. 






m 


v m ,o 


5 2 v 


8*v 


8+v 


5°v 


-1 

1 
2 
3 
4 
5 


0. 9166802 
. 8000000 
. 9203245 

1. 2164087 

1. 5985781 

2. 0228433 
2. 4711277 


. 2370047 
. 1757597 
. 0860852 
. 0420958 
. 0240192 


-. 0612450 
-. 0896745 
-. 0439894 
-. 0180766 


-. 0284295 
+ . 0456851 
+ . 0259128 


+ . 0741146 
-. 0197723 



The fourth differences are not numerically smaller 
than the third differences, and the first entry in the 
fifth-difference column is very much larger numeri- 
cally than the second one in the same column. Such 
a pattern is a warning that the interval is too large. 

<d) The Effect of Using a Higher-Order Approximation Formula 

Let us compare the results in examples 1 and 4. 
The coding for example 4 is somewhat more compli- 
cated than that for example 1, but since only 1,472 
lattice points were used in example 4 compared with 



2,000 in example 1, the over-all amount of work is 
about the same in both cases. The result after 100 
steps shows that the higher approximation formula 
gives very much better results. To secure a maxi- 
mum error of 0.00008 in v m>n with the approximation 
used in example 1, it would have been necessary to 
use an interval h of about 0.009; hence 170 times the 
amount of work would have been necessary. How- 
ever, if the p 2 correction were applied, the compara- 
tive results would not be quite so unfavorable to the 
simpler approximation. Assuming that a p 2 correc- 
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tion were applied to two computations by the simpler 
approximation, and a p 4 correction to results of the 
modified fourth-order approximation, the latter 
would still give significantly better results — one 
additional decimal place, in fact. To secure com- 
parable accuracy by the simpler approximation, it 
would be necessary to multiply the amount of work 
by the factor 7 10. The conclusion is inescapable 
that the higher approximation is worth while, in 
cases where the boundary conditions do not introduce 
singularities in the higher derivatives, and when the 
coding problem is manageable. 



The author acknowledges gratefully the many 
constructive suggestions which were given by Dr. 
Fritz John during the progress of the study. 
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7. Appendix 

Definition. We define z (m,M;n,N) to be a region in the 
z,£-plaiie, which contains the region 

(?n - M) h <x < (m + M) h; nk<t<(n+ N) k. 

7.1. Expressions for Derivatives 8 in Terms of 
Differences 



dtu 1 
dx' 



4[( 52 -T2 64 +i 56 ) w -]- 



560 dz 8 ' 



(37) 



dx 



^4[(H s °K»>'[ 



7 d*u x „ h 2 d«>tt- „- 



240 dz 8 4200 dx 10 



(38) 



d 6 u„ 



dx* 



n -=±8«u m . n -h*[ 



l<**U Xi , n 7h 2 d™U„ n hWu„ n - 



4 dx 8 720 dx 10 ' 184805a- 12 



(39) 



where the points (x if nk) are in Z(w,3;n,0), j=l,2,. . .,6. 



7 This is on the assumption that after the p 2 correction, the leading term of 
<Tm. n has the factor /i 3 . Since the work is inversely proportional to h* for the same 
X, improving the entries by a factor of At means doing M times the amount of 
work. 

8 See chapter VII, reference [3]. 



If h is sufficiently small, successive terms of (37), (38), and 
(39) are of a lower order of magnitude than preceding ones, 
and if terms involving h in the numerator are dropped, then 
the truncation error is of the order of magnitude of the first 
term neglected. 

7.2. Relations Between Derivatives in the /- and 
./-Directions 



By differentiating (I) we can establish: 

d 2 U== d 2 /du\ df = d*u d 2 f df 
dt 2 ~ dx 2 \dt)^~ A * ~ *~* ~*~ ^ 2 



dt 0.r 4 dx 2 



dt 



d 2 u = d*u df 
dxdt dx 3 dx 

&u _d% d*f d 2 f 

?w5+ /31 > ^ 31 dx^dxdt 



dxdt 2 dx 
d 3 u d & u 



dt 3 dx 6 



P32 ; P&. 



d'f cPf d*f 
~dx 4 ^~dP^~dldx* 



(40> 
(41> 
(42) 
(43) 



Iii the above, Pu and P 32 are functions of a; and t. Where not 
otherwise specified, they will be understood to correspond to- 
x — mh, t = nk. If f(x,t,u) involves u, the derivatives of/ may 
involve various derivatives of u with respect to x. It will be 
convenient in this section to change the notation introduced 
previously, and to define f x = df/dx, with t,u held fixed; f t — 
df/dt, wifh.r,// held fixed; /u=d//dw, with x,t held fixed witn sim- 
ilar definitions for / z , «,/*,«, etc. These derivatives will usually 
be required at x = and t = nk; hence evaluation of the deriva- 
tives at (0,nk) will be implied unless otherwise specified, or 
readily understood from the context of the analysis. 



7.3. The Boundary Conditions at x = 
Examples in Section 5 



for the 



Since u- m ,o has not been defined by the differential system, 
the expression (47) must be modified when m = 0, or what is* 
equivalent to it, u-i f „ must be defined. Similarly u-\ >n and 
u-2, n must be defined for (5) and (32). Consider the condi- 
tion (dwo,»/daO=0. It implies 



dt p ~ l dx 



Now from (41) and (44) 



=0, p>2. 



+/*■ 



(44) 





_d 2 M O ,n_C) 3 W0,n 

dtdx dx 3 


Hence 


d*u , n . 
dx 3 Jx ' 


Similarly, 


from (6.23) and (6.30) 


where 


d Z U , n d 5 U ,n 

dt 2 dx dx 5 



(45> 



+ P31, 



P =r*f+#r\ . 

3,1 Ld&^dxdtJo.t 

If / involves u, P31 will involve d 2 u/dx 2 . An explicit 1 ex- 
pression for P31 in terms of partial derivatives will therefore 
be useful. It is given below. 



-/xxx + 4/ x , m t— ^-+ f'fx,u~~fz'fu~\~fx,i 



--*fz.u 



d 2 u Q , 
dx 2 



^•--sr 2 ' (46> 
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where 



O n —fxxx +/'/x, u — fx'fu +/*,!• 



(47) 



Using (45) the Maclaurin series around (0, nk) yields 



U m ,n = UQ t n + 



m 2 h 2 d 2 tt ,n w 3 ^ 3 



dx 2 



6 p=4 p! 



m p /i p d p u ,n 



dx v 



(48) 



Ifwesetra=l, 2 and use (46), we can obtain d 2p u 0jn /dx 2v to 
an accuracy comparable to approximations used over the 
rest of the region. The same results are obtained if we use 
(48), to solve for u~\, n and u- 2 , n (if the latter is needed). 
This artificial extension of the region to negative values of 
m is convenient for numerical treatment, and is justifiable 
wherever (48) exists and the required derivatives are bounded. 
Thus we can write 



m z h z . m 5 h 5 d 5 Uo t . 



-f. 



60 dx 5 



(49) 



If h is sufficiently small, and terms involving powers of 
h 5+v ,p>0 are neglected in (49), the truncation error for 
u- m ,n is of order h 5 . Thus T 2t o f n, defined in (7a), has a 
term in h 3 . 



The third term on the right-hand side of (49) was also 
dropped in example 4, since its magnitude would have 
affected only the fifth decimal place at any point. This 
term could have been obtained, if desired, by using (46) 
and (37), and then setting m=l and m = 2 in (49), to obtain 
three linear equations for the unknowns u-i,u- 2 , and h 5 d 5 u/dx 5 . 

7.4. The Boundary Condition at the Terminal 
Points (x a , t) 

The difference equation defined in (4) needs no special 
treatment for the boundary conditions at x a , where a = sm, 
since (s — l,n) is the last lattice point at which v m ,n is gen- 
erated, and 8 2 v s -i fn is fully defined. However, when the 
fourth order approximation is used, as in (5), an expression 
is required for v s+ltn , or for 8 2 u s>n . Since u(x a ,t) = B(t) and 
all its derivatives in the ^-direction are assumed to be known, 
the differential eq (1) can be used to obtain the second 
partial derivative in the ^-direction, in terms of du/dt and 
f(x,u). The relations between derivatives and differences can 
then be used to obtain 8 2 u Sfn . In examples 4 and 5, the 
fourth difference in the x-direction was essentially zero at 
x a ; hence the last fourth difference was replaced by 5 4 w s _i >n . 

Los Angeles, California, September 18, 1951. 
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